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We present an efficient and general method to compute vortex-pinning interactions - as appear in 
neutron stars, superconductors, and trapped cold atoms - at arbitrary separations using real-time 
dynamics. This method overcomes uncertainties (associated with matter redistribution by the 
vortex position and the related choice of ensemble) that plague the typical approach of comparing 
energy differences between stationary pinned and unpinned configurations: uncertainties that prevent 
agreement in the literature on the sign and magnitude of the vortex-nucleus interaction. We 
demonstrate and validate the method with Gross-Pitaevskii-like equations for the unitary Fermi gas, 
and demonstrate how the technique of adiabatic state preparation with time-dependent simulation 
can be used to calculate vortex-pinning interactions in fermionic systems such as the vortex-nucleus 
interaction in the crust of neutron stars. 



PACS numbers: 97.60.Jd 26.60.-c 26.60.Gj 21.60.-n 

Vortex-pinning interactions play an important role in the 
dynamics of various condensed superfluid systems such 
as superconductors [T], trapped cold-atom gases and 
possibly neutron stars f3], where the angular momentum 
carried by vortices can have an observable impact. Pulsar 
glitches, for example - sudden increases in the rotation 
frequencies of neutron stars - are theorized [4 to arise 
from a sudden macroscopic unpinning of vortices: In equi- 
librium, the superfluid and non-superfluid components of 
a pulsar rotate at the same angular frequency. The pulsar 
loses angular momentum through magnetic radiation, and 
the crust slows down gradually, reducing the pulsation 
rate. To maintain equilibrium, the superfluid must also 
release angular momentum by diluting the vortex concen- 
tration, but the presence of pinning sites (nuclei, lattices 
sites, defects, etc.) may arrest the vortex motion; stress 
would build until a large number of vortices rapidly unpin, 
dilute, and transfer their angular momentum to the crust, 
rapidly increasing in the pulsation rate - the glitch. 

Despite almost 40 years, the feasibility of this mecha- 
nism is still poorly understood. The conventional picture 
has the angular momentum stored by the neutron su- 
perfluid in the crust (which is well approximated by the 
same unitary Fermi gas (UFG) produced in cold-atom 
experiments [5]), with pinning provided by nuclei held 
in a lattice by the electrostatic (Coloumb) interaction. 
Recent results suggest that the crustal neutrons may not 
support enough angular momentum to explain observed 
pulsar glitches [6_ , in which case pinning in quark matter 
phases in the core may play a role [7]. In either case, a 
reliable technique for calculating vortex-pinning interac- 
tions is key. In this paper, we present a dynamical method 
for determining the sign and strength of vortex-pinning 
forces in superfluids, and demonstrate that this method 
can be directly applied to unambiguously calculate the 
vortex-nucleus interaction using time-dependent density 



functional theory (TDDFT) for nuclear matter. 

Due to the importance of pinning on glitch phenomenol- 
ogy, several attempts have been made to calculate the pin- 
ning force in nuclear matter from underlying microscopic 
models. The earliest calculations ^ used a Ginzburg- 
Landau (GL) calculation to estimate the pinning energy 
as a function of density. In [51 the GL framework was used 
to give a detailed picture of the pinning/unpinning pro- 
cess: their calculation includes an estimate of the energy 
as a function of displacement allowing for an estimate of 
the pinning force. The next advance was the use of a local 
density approximation [TU] with Wigner-Seitz cells and 
Gogny [XT' and Argonne fT^ interactions, which gives 
similar density-dependant pattern of pinning/unpinning 
as [9], but smaller by almost an order of magnitude. 

Unlike vortices in weakly coupled BCS superfluids, vor- 
tices in dilute neutron matter (and the UFG) displace a 
substantial amount of matter from their core |13| . There- 
fore, comparing "energies" of stationary configurations 
with a nucleus on the core of a vortex and a nucleus away 
from the vortex is confounded by a choice of ensemble: 
should one fix the number of neutrons or the chemical 
potential in a finite simulation volume? 

Computing stationary configurations is also computa- 
tionally expensive - especially given the high degree of 
precision required to render meaningful energy differences 
- and simulations to date have required a high degree of 
symmetry. For example, recent self-consistent calcula- 
tions [T3HTB] extract the binding energy of a vortex on 
a single nucleus in a cylindrical geometry, which is not 
the same as the pinning energy. In particular, the con- 
clusion of |15] that the pinning force is repulsive (glitches 
would thereby requiring interstitial pinning) was ques- 
tioned by [T7] , and a different set of calculations using the 
local density approximation suggests that pinning is at- 
tractive over a substantial region in the inner crust |14([T5]. 
Moreover, nearby vortices and the Casimir effect can sig- 
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FIG. 1. (Colour online) Deflection of a vortex in the ETF model of trapped dilute neutron matter as a UFG by a repulsive 
(left panel) /attractive (right panel) pinning potential Vpin(r) = ±3.5MeV/[l + exp(r/fm — 7.5)] moving on a straight line from 
left to right at a constant subsonic velocity v ~ O.lCs. The trajectory of the vortex is shown by the (black) curve and the 
relative separation vector between the pinning site and the vortex core is shown as thin (white) lines for select times connecting 
the corresponding dots on the trajectories. Initially the potential displaces the bulk superfluid, carrying the vortex to the 
right/left. Once the potential overlaps with the vortex, the vortex rapidly moves down/up - (almost) perpendicular to the 
force. In the frame shown on the left, the pinning site is just to the left of the centre {x ~ — 2.5fm) and the vortex is moving 
(almost) perpendicular along the edge of the pinning potential. After the potential has passed through, the vortex orbits in 
a counter-clockwise circle direction due to boundary effects from the trap that can be quantitatively described in this sharp, 
flat trap by placing an image vortex outside of the potential to cancel the tangential current at the boundary: this induces a 
counter-clockwise super-flow Vs in Eq. |[T]). The geometry of the right simulation is such that the potential carries the vortex 
around almost the entire trap: This extended interaction allows the pinning potential to excite phonons in the system visible as 
ripples in the circular trajectory. 



nificantly polarize a nucleus - an effect absent in simple 
cylindrical geometries - dramatically changing the nature 
of the nuclear pinning sites and disrupting the regularity 
of the nuclear lattice [12] • 

Characterizing the nuclear-pinning interaction will thus 
require self-consistent calculations in 3D without symme- 
try constraints: currently not feasible with techniques 
that require highly-accurate stationary states, but feasi- 
ble with TDDFT where the algorithms are substantially 
more efficient and exhibit better scaling properties. The 
application of TDDFT for studying fully 3D superfluid 
dynamic has been demonstrated in [20_ and scales well 
to massively parallel supercomputers. We now present a 
qualitatively new approach for calculating vortex-pinning 
interactions, unencumbered by the aforementioned issues, 
and utilizing only real-time dynamics. 

The idea, similar to the Stern-Gerlach experiment, is 
to observe how a vortex moves when approached by a 
nucleus. To zeroth order, the sign of the interaction is 
determined qualitatively by the direction of the motion 
(Fig.[l]); with a more careful inspection, one can extract 
the force-separation relationship F{r) (Fig. [2]). 

We validate our procedure using a dynamical extended 
Thomas- Fermi (ETF) model ^21, - ,24) . equivalent to a 



Gross-Pitaevskii equation (GPE) for bosonic "dimers" 
niB = 2m of fermionic pairs, with an equation of state 
tuned to consistently fit both quantum Monte Carlo 
(QMC) and experimental results [53]. Despite the compu- 
tational simplicity of the ETF model, it has been demon- 
strated to quantitatively reproduce low-energy dynamics 
of both UFG experiments [22] and fermionic density func- 
tional theory (DFT) simulations p4|; in the latter case, 
the ETF also characterizes numerical effects such as fi- 
nite size and discretization artifacts. The UFG should 
also qualitatively model the dilute neutron superfluid in 
the crust of neutron stars |25| due to the large neutron- 
neutron scattering length a„„ — 18.9fm [2^. Thus, by 
using a physically motivated model of the nuclear pairing 
potential [IT,, we anticipate these ETF calculations to 
faithfully simulate future fermionic TDDFT simulations. 

To gain some intuition for the vortex-nucleus interac- 
tion, consider the phenomenological Hall-Vinen-Iordanskii 
(HVI) equation (see [27| for a discussion) for a vortex in 
2D: 

Mfv - fqp = Psii X in, - Vs) + Fy. (1) 

Here is the position of the vortex, the force Fy is per 
unit length along the vortex, ps is the number density of 
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the "background" superfluid, k = 2TThz is the quantized 
vortex circulation, and Vs is the "background" superfluid 
velocity. This equation should only be taken as an intu- 
itive guide since terms on the left-hand side are ill-defined. 
The "mass of the vortex" M, for example, depends strongly 
on the way it is measured |5S] , and the force fqp due to 
excited phonons has significant memory effects. 

One can argue that the left side of this equation is of 
higher order than the right side for small velocities of 
the pinning site, and, as we shall verify in Fig. [2| can 
be neglected for the application of vortices in neutron 
stars where the pinning force is quite weak and moves 
slowly. This leaves the well-established Magnus relation- 
ship psK X (r\, — iT,) « —Fy relating the force Fy applied 
to the vortex and its perpendicular velocity relative to 
the background superfluid velocity Vg. Thus, by observ- 
ing the dynamical deflection of a vortex from a nuclear 
pinning site, one can directly extract the direction and ap- 
proximate magnitude of the vortex-nucleus force without 
requiring a subtle subtraction of energies. 

In small systems, the Magnus relation can only be used 
to estimate the magnitude of the force since the superfluid 
density ps and velocity Vs are not precisely defined, though 
reasonable estimates can be obtained. With an external 
pinning potential Vpia{fpin — , however, one can directly 
and unambiguously calculate the force on the pinning 
site: 

f pin = - / d^X ^^P'-^y-'^^^ plx). (2) 

In the nuclear context where neutrons are present in 
the both the pinning site (the nucleus) and the superfiuid 
medium, the force can be obtained in two ways: 1) Eq. ^ 
can be directly applied to a Coulomb potential that traps 
the protons - this will be the force that the vortex ex- 
erts on the nuclear lattice - or, 2) one can estimate the 
force using Newton's law -Fpin = rripi^api^ for a dynamic 
pinning site comprising protons and entrained neutrons. 
The position of the pinning site can be unambiguously 
defined as the center of mass (CM) of the protons, and 
the effective mass rripin can be estimated limiting forms 
in homogeneous matter |29j. 

What remains is to prepare the initial conditions 
with a vortex and nucleus interacting at various dis- 
tances. We can do this using adiabatic state prepara- 
tion [30 , thereby eschewing the need for diagonalizing 
large matrices to find stationary initial conditions. The 
idea is to adiabatically evolve in real-time a state of 
some solvable system to a desired initial state in the 
system of interest. For example, starting with a non- 
interacting harmonically trapped gas, we can evolve ei- 
ther the ground state 'Fgs exp(— TOBa;r^/2/i), or an 
exact "Landau level" 'Fs oc {x + iy — S) exp(—mBi^r'^ /^h) 
(stationary in a rotating frame) with angular momentum 
Iz = Nh/{\ + rabjS^ jK). We adiabatically evolve the sys- 
tem while turning on the interactions and adjusting 



the trapping potential Vt = (1 — s)Vfro -I- sVtrap with 

a C°° function s{t/T) = / that goes from to 1 

over a characteristic time T chosen to be longer than any 
intrinsic timescale in the system: 
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(The parameter a controls the smoothness and we find 
empirically that values 1 < a ^ 6: ^^/"^ work well.) 
From '^GS we can generate the ground state, and from 
!f'5=o we can generate a single vortex in the centre of the 
trap, both to high precision. Producing a state with a 
single orbiting vortex off-axis vortex proves to be more 
challenging - one often finds incommensurate "phonon 
sloshing" in addition to an orbiting vortex (similar to the 
final state in the right of Fig. [ij. 

To accurately measure the vortex-pinning interaction 
one can choose as a final potential Vt = T4rap + ^pin: 
an axially symmetric trap of suitably flat bottom and a 
pinning potential in the center. By generating a configura- 
tion with a vortex orbiting in a circle at radius r, we can 
use Eq. ([2| to calculate the force exerted on the centrally 
located pinning potential: axial symmetry ensures that 
this is precisely the vortex-pinning force at separation 
r. We use this procedure within the ETF model to ac- 
curately calculate the "stationary" (in a rotating frame) 
vortex-pinning interaction shown in Fig. [2] 

This procedure is not ideal with the adiabatic state 
preparation due to the sloshing. Instead, one can prepare 
a central vortex tf'5=0j ^■iid then directly probe the force 
using a "dynamic" simulation as shown in Fig.[l]for a large 
range of interaction. (The perpendicular motion of the 
vortex makes it somewhat difficult to probe very short- 
range interactions with this approach.) The feasibility 
of this approach is verified in Fig. |2] where we directly 
compare it with the vortex-pinning interaction extracted 
by the stationary method. 

The present demonstration has been limited to quasi- 
2D simulations. The procedure will work just as well in 
fully 3D simulations. New effects such as the bending 
of a vortex line when approached by a pinning site can 
just as easily be analyzed: the vortex line will either 
be repelled by the pinning site - bowing out to avoid 
it - or will be sucked in. Again, qualitative dynamical 
behaviour coupled with the techniques demonstrated here 
can be used to extract the 3D characteristics of the vortex- 
pinning interaction. We have considered here only moving 
the pinning site, but one could also consider manipulating 
parts of the vortex with pinning potentials, dragging the 
pinned vortex along a trajectory instead. In simulations 
with realistic nuclei, the vortex-nucleus interaction will 
also excite and deform the nucleus - significantly affecting 
the vortex-nucleus interaction. It is conceivable also that 
a vortex lines could break and attach to various nuclear 
defects like rods of plates: the dynamics of such broken 
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FIG. 2. (color online) Here we demonstrate the consistency of 
dynamically extracting the vortex-pinning force for the nuclear 
pairing potential tllj V"pin(r) = 0.75 MeV/[l + exp(r/fm- 7.5)] 
at densities p ~ 0.045 fm""^ k, 0.28psat. The triangular (blue) 
points come from the computationally expensive "stationary" 
method, while the solid (green) curve come from using the 
"dynamic" real-time evolution analogous to that shown on the 
left panel of Fig. [l] The dotted (red) curve shows the Magnus 
estimate for the force ([T]) using a Thomas- Fermi approximation 
for ps and estimating Vs from the image vortex [31| . The 
double curves represent the interaction as measured from the 
incoming and outgoing pinning site. 



vortex lines may also play a important part in explaining 
neutron star glitches. 

We close with a brief analysis of the time-evolution. 
The complex scalar field If" obeys an evolution equation 
of the form 



IrriB 



(4) 



where yeffl!^] is an effective interaction that depends non- 
linearly on and on the external trapping and pinning 
potentials. Consider the quasi-two-dimensional problem 
where coordinates may be expressed as complex numbers 
z = X + iy: A singly-wound vortex at location may be 
described by the field tf'(z) = [z — Zy)f{z), where f{z) is 
a smooth complex-valued function that we assume has 
no roots in the immediate vicinity of the vortex. The 
evolution equations may then be expressed as follows 



ihf 



2m B 



_ ihzyf - h^{dxf + \dyf)/mB 



Z Zd 



(5) 



The left side is smooth, hence the pole on the right side 
must cancel with a root in the numerator, giving us an 



explicit expression for the vortex velocity 

^ -i9x/ + dyf 



Zv = \rv\x + iK 



■niBf 



(6) 



This expression can be written as the "Magnus" relation 



K X (fy — V) 



2mBP 

z—z^ 

where the local "superfluid velocity" v is 



(7) 



(8) 



Note that the relationship ([t]) is exact. The meaning of 
Ps in the HVI equations is not clarified since p cancels 
in ^ . Unfortunately, although v is precisely defined and 
corresponds to Vs in some situations, one cannot generally 
make the correspondence Vg = v. In particular, doing so 
yields results that differ by as much as 50% in Fig. [2] 

The quantities appearing in Eq. ([T]) are related to long- 
range momentum transfers and boundary effects, and one 
must thus be content with reasonable estimates for ps 
and Vg , for example, from the average behaviours of the 
relevant quantities near but outside of the vortex core. 
As Fig. |2] demonstrates, however, the Magnus relation 
is suitable for extracting the sign and magnitude of the 
interaction. This also provides an explicit check that the 
force evaluated using our procedure is what appears on 
the right side in Eq. ([T| governing the vortex dynamics. 

Conclusion: We have described how to use time- 
dependent density functional theory (TDDFT) to effi- 
ciently and unambiguously calculate vortex-pinning in- 
teractions from real-time dynamical simulations of super- 
fluid systems. We have demonstrated with an extended 
Thomas-Fermi (ETF) model of the unitary Fermi gas 
(UFG) that this approach can be applied to calculate 
the vortex-nucleus interaction using nuclear TDDFTs to 
model the crust of neutron stars. While we considered 
only quasi-2D systems here, the size of the problem, the 
magnitude and accuracy of the force extraction, and the 
use of pure real-time dynamics ensure that full 3D simu- 
lations of realistic fermionic TDDFTs will be within the 
capabilities of next generation supercomputers. A resolu- 
tion to the puzzle of pulsar glitches will requiring more 
than just extracting the vortex-nucleus interaction, but 
with this real-time method, this crucial step will soon be 
within reach. 

These techniques can also be directly applied to systems 
of trapped ultra-cold atoms in a variety of geometries, for 
example, to explore vortex pinning on optical lattices. In 
particular, the close approximation of the neutron super- 
fluid by the UFG suggests that cold atom experiments 
might also be able to shed light on the glitching puzzle. 
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